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We develop and calibrate a characteristic waveform extraction tool whose major improvements 
and corrections of prior versions allow satisfaction of the accuracy standards required for advanced 
LIGO data analysis. The extraction tool uses a characteristic evolution code to propagate numerical 
data on an inner worldtube supplied by a 3+1 Cauchy evolution to obtain the gravitational waveform 
at null infinity. With the new extraction tool, high accuracy and convergence of the numerical error 
can be demonstrated for an inspiral and merger of mass M binary black holes even for an extraction 
, worldtube radius as small as R = 20M. The tool provides a means for unambiguous comparison 

qq ■ between waveforms generated by evolution codes based upon different formulations of the Einstein 

equations and based upon different numerical approximations. 

PACS numbers: 04.20Ex, 04.25Dm, 04.25Nx, 04.70Bw 
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I. INTRODUCTION 

The strong emission of gravitational waves from the inspiral and merger of binary black holes has been a dominant 
motivation for the construction of gravitational wave observatories. The computation of the precise details of the 
waveform by means of numerical simulation is a key theoretical tool to enhance detection and allow useful scientific 
interpretation of the gravitational signal. See 1] for a review of the accuracy required of numerically generated 
waveforms to fully complement the sensitivity of the LIGO and Virgo @ observatories. However, the waveforms 
are not easy to extract accurately. The radiation falls off as 1/r so that although it asymptotically dominates near 
field gravitational effects it is nevertheless small and can be contaminated by numerical error. It is common practice 
for the Cauchy codes used in simulating the binary black hole problem to introduce a large but finite artificial outer 
boundary. A combination of linearized and far field approximations are then used to extract the waveform from data 
on a smaller worldtube, which ideally is causally isolated from the outer boundary. Such perturbative wave extraction 
at a finite distance, rather than at null infinity which more faithfully represents the idealization of a distant antenna, 
introduces systematic errors associated with the effects of gauge, nonlinearities, nonradiative near fields and back 
reflection. (See 0, HI for analyses of waveform errors arising from perturbative extraction at a finite distance.) An 
alternative approach called Cauchy- characteristic extraction (CCE) provides a fully nonlinear interface between 
Cauchy and characteristic codes which utilizes the characteristic evolution to extend the simulation to null infinity, 
where the waveform is computed. An earlier implementation of CCE has recently been applied to extract waveforms 
from binary black hole simulations p| Q , from rotating stellar core collapse [l(| and to explore the memory effect [ll[ . 
In this work, we present details and tests of a redesigned CCE module whose accuracy and efficiency has undergone 
major improvement. The module has been designed to provide a standardized waveform extraction tool for the 
numerical relativity community which will allow CCE to be readily applied to a generic Cauchy code. 

The first attempts to simulate collisions of black holes by Hahn and Lindquist [12| . and then by Smarr et al 13], 
were hampered by both a lack of computing power and a proper understanding of the mathematical formulation of 
Einstein's equations required for a stable numerical solution. Their work formed the impetus for the Binary Black 
Hole Grand Challenge, which was formed to take advantage of the increasingly powerful computers introduced in the 
1980's. The main results of the Grand Challenge were limited to the axisymmetric head-on collision of black holes 
and the gravitational collapse of rotating matter However, the standard Arnowitt-Deser-Misner [lj| formulation 
of the Einstein equations adopted by the Grand Challenge had instabilities at the analytic level which limited more 
general binary black hole simulations to the premerger stage. Only with new formulations was a full inspiral and 
merger successful, first by Pretorius [l6j using the harmonic formulation, and soon after by Campanelli et al (l7j and 
Baker et al [HI using the Baumgarte-Shapiro-Shibabta-Nakamura formulation [lj|[20[. Numerous groups now have 
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codes which can simulate this binary inspiral problem by evolving the Cauchy problem for Einstein's equations. 

In CCE, the Cauchy evolution is used to supply boundary data on a timelike inner worldtube to carry out a 
characteristic evolution extending to future null infinity where the waveform can be unambiguously computed 
using the geometric methods developed by Bondi et al [2l|, Sachs [22j and Penrose [23j |. This initial-boundary value 
problem based upon a timelike worldtube [241 ] has been implemented as a characteristic evolution code, the PITT null 
code [25L l26j. which incorporates a Penrose compactification of the space-time. It computes the Bondi news function 
at which is an invariantly defined complex radiation amplitude N = N® + iN®, whose real and imaginary parts 
correspond to the time derivatives and dth® of the "plus" and "cross" polarization modes of the strain incident 
on a gravitational wave antenna. The error in the PITT code was tested to be second order convergent in analytic 
testbeds ranging from the perturbative regime [13] to highly nonlinear single black hole spacetimes [26j]. One of the 
successes of the Grand Challenge was the successful application of the code to generic single black hole dynamical 
spacetimes [28M3TI ] . For a review, see [32| . 

The propagation of gravitational waves to from an astrophysically realistic source using the PITT code has in 
the past been limited to the simulation of an imploding neutron star using a fluid dynamic code incorporated into 
the characteristic code [33|, H3] • These simulations were restricted to the axisymmetric case because of computational 
demands arising at the center of the star. For such systems, CCE offers a way to combine the strengths of the Cauchy 
and characteristic approaches. Recently this combined approach has been applied to extract the waveform from the 
fully 3-dimensional collapse of a rotating star [lOj . A global characteristic simulation of the full inspiral and merger 
of a relativistic binary system is not possible because of the interior caustics formed by gravitational lensing. But the 
application of CCE to this system has been shown to be possible [§|, H|. 

The error in CCE arises from three independent sources: (I) the Cauchy evolution; (II) the worldtube module; 
and (III) the characteristic evolution to J'^ and the computation of the waveform. 

(I) Errors in the Cauchy evolution can arise from numerical approximations, improper boundary treatments, 
extraneous radiation content in the initial data, instabilities and bugs. Errors introduced at the outer grid boundary 
present a special problem for BSSN formulations for which there is no theoretical understanding of the proper boundary 
condition. The standard practice is to extract the waveform at a finite worldtube which is large enough to justify a 
far field approximation but which is causally isolated from the outer boundary during the simulation. For example, 
perturbative extraction at r = lOOAf would require that the outer boundary be at r > 500M for a t w 400M 
simulation. We have designed the new CCE module so that it can be applied to a generic Cauchy code with extraction 
radius as small as r — 20M. However, since any universally applicable extraction module must be designed to be 
independent of error introduced by the Cauchy code, the extracted waveform cannot be any more reliable than the 
Cauchy code. 

(II) The main improvement described and tested in this paper is a complete overhaul of the worldtube module, 
which converts the output of the Cauchy evolution to boundary data on an inner worldtube for the characteristic 
evolution. The prior version of this module, which was used in the first applications of CCE to obtain binary black 
hole waveforms [8l4lOl]. contained inconsistencies and bugs which prevented clean convergence tests. We have corrected 
this worldtube module so that the present version exhibits clean convergence to which Richardson extrapolation can 
be applied to produce waveforms whose numerical error due to CCE is extremely small. In addition to improvement 
in consistency and accuracy, we have also redesigned the module to be more efficient and user friendly. These revisions 
are described in Appendix |A"1 

(IIII) In addition to thoroughly scrutinizing the PITT null code for bugs, we made several major modifications. In 
previous applications requiring very high resolution, such as the inspiral of a particle into a black hole [35| , there was 
excessive short wavelength noise which affected the quality of the simulation. In addition, in [1,0 it was reported that 
one of the equations governing calculation of the waveform at had to be linearized in order to obtain reasonable 
behavior. These problems have been eliminated as a result of the modifications described in Appendix [Al 

In See's. [IT] and IIII1 we review the formalism underlying characteristic evolution and the computational structure 
of the PITT code. We include enough details to make clear the difficulties underlying extraction of an accurate 
waveform at ^ + and to explain the code modifications that have been made. We also demonstrate how the use of 4th 
order accurate angular derivatives improves the previous test results of CCE presented in [5]. In Sec. IIV1 we describe 
the design of the new worldtube module, how it treats the Cauchy-characteristic interface and how it can be readily 
applied to a Cauchy evolution. 

In Sec. El we test the new extraction tool on the Cauchy evolution of the inspiral and merger of two equal-mass, 
non-spinning black holes. We show that CCE can now be carried out for a worldtube radius as small as 20M for 
a mass M binary system, for which perturbative extraction would not be meaningful, and which was not possible 
with the prior implementation of CCE. Convergence tests now demonstrate clean second order global accuracy of the 
evolution variables. 

The waveforms are only first order accurate as a result of the asymptotic limits required at .y + . However, the clean 
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first order convergence of the waveform now allows application of Richardson extrapolation to obtain higher order 
accuracy. In this way, in Sec. IVI[ we construct a third order accurate waveform, which was not possible with earlier 
versions of CCE. 

The ability to apply Richardson extrapolation to CCE waveforms makes it possible to show that their numerical 
error satisfies the standards required for application to advanced LIGO data analysis. The first derivation [36[ of the 
accuracy required for numerically generated black hole waveforms to be useful as templates for LIGO was carried out 
in the frequency domain, in which the interferometer noise spectrum is calibrated. There are two separate criteria: 
one ensures that the error in the model waveform does not impact wave detection and the other ensures that the 
error does not impact the scientific content of the signal. These criteria both depend upon the noise spectrum of the 
detector in a way not easily applied to a numerical simulation. This has recently prompted a translation of these 
requirements into the time domain in which the waveforms are computed P, l37l |38| , so that they can be readily 
enforced in practice. In Sec. I VIII we show that the numerical error introduced by CCE satisfies these time domain 
criteria for an advanced LIGO detector. We also analyze the error introduced by the choice of initial data, which has 
a dependency upon the size of the extraction worldtube. 



II. CHARACTERISTIC FORMALISM 



The characteristic formalism is based upon a family of outgoing null hypersurfaces emanating from an inner world- 
tube and extending to infinity where they foliate ^ + into spherical slices. We let u label these hypersurfaces, x A 
(A = 2, 3) be angular coordinates which label the null rays and r be a surface area coordinate. In the resulting 
i° = (u,r, x A ) coordinates, the metric takes the Bondi-Sachs form [2~il I22I] 

ds 2 = -( e 2p — - r 2 h AB U A U B ) du 2 - 2e 20 dudr - 2r 2 h AB U B dudx A 



r 

+ r 2 h AB dx A dx B , (2.1) 

where h AB h B c = and det(h AB ) = det(q AB ), with q AB a unit sphere metric. In analyzing the Einstein equations, 
we also use the intermediate variable 

Q A = r 2 e~ 2f3 h AB U B r . (2.2) 

Because the Bondi variable V — 0(r 2 ) at the code is written in terms or the renormalized variable W = (V—r)/r 2 . 
Here W = for the Minkowski metric in null spherical coordinates. 

The PITT null code employs a spherical grid based upon an auxiliary unit sphere metric q AB , with associated 
complex dyad q A satisfying q AB — ^ ( K q A q B +q A q B ). The Bondi-Sachs metric h AB induced on the spherical cross- 
sections can then be represented by its dyad component J = h AB q A q B /2, with the spherically symmetric case 
characterized by J = 0. The fully nonlinear h AB is uniquely determined by </, which is the principle evolution variable. 
The determinant condition implies that the dyad component K = h AB q A q B /2 is determined by 1 = K 2 — J J. We 
also introduce spin- weighted fields U — U A q A and Q = QaQ A , as well as the complex spin- weight operators 8 and 
9 [39| which represent the angular derivatives. Refer to (ioj for details regarding numerical implementation. 

In this formalism, the Einstein equations decompose into hypersurface equations, evolution equations and conser- 
vation conditions on the inner worldtube. As described in more detail in @, HH, the hypersurface equations take the 
form 

P,t = N P \J\, (2.3) 
{r 2 Q), r = -r 2 (dJ + dK), r + 2r 4 d(r- 2 f3) r + N Q {J : f3], (2.4) 

U r = r- 2 e 2 PQ + Nu[J,p,Q], (2.5) 

V r = ^-e^K-ePmeP + ^-r~ 2 (r 4 (W + 3U)) r + N W [J, /3,Q,U], (2.6) 

where 

n = 2K- mx + i(8 2 j + g 2 J) + -^(SJcV - B.mJ) (2.7) 

is the curvature scalar of the 2-metric h AB . The evolution equation for J takes the form 

2(rJ) ur -(r~ 1 V(rJ) r ) r = 

-r- 1 {r 2 W) r + 2r-VciV - (r~ l V) p J + Nj[J, J„,/3, Q, U, W], (2.8) 
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where Np[J], Nq[J,/3], Njj[J, f3,Q], Nw[J,@,Q,U] and Nj[J,J tU ,ft,Q,U,W] are nonlinear terms which vanish for 
spherical symmetry and can be constructed from the hypersurface values of the variables appearing in their argument. 
Expressions for these nonlinear terms as complex spin- weighted fields and a discussion of the conservation conditions 
are given in Q. The hypersurface equations have a hierarchical structure in the order [J, /3,Q,U,W] such that the 
right hand sides, e..g. Np[J] only depend upon previous variables and their derivatives intrinsic to the hypersurface. 
The finite difference grid used in the code is based upon the compactified radial coordinate 

(2.9) 



Re + r 

so that x = 1 at J f+ . Here Re is a parameter which in the CCE module is chosen as the radius of the extraction 
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worldtube as determined by R 2 = 8ijX l x^ in terms of the Cartesian coordinates x % used in the Cauchy evolution code. 



The auxiliary variables 

u = dJ, 6 = 5/3, k = dK (2.10) 

are also introduced to eliminate all second angular derivatives. In certain applications this has been found to give rise 
to increased accuracy by suppressing short wavelength error |4lj . 

The finite difference scheme for integrating the hypersurface and evolution equations has been described in (26l . l4ll . 
|42| . Except for the start-up procedure described in Sec. IIV1 we follow this scheme with two modifications. First, the 
finite difference approximation for the 9-operators is increased from 2nd order to 4th order accuracy. This can be 
expected to give better angular resolution but does not affect the overall 2nd order accuracy implied by the radial 
and time integration schemes. Second, when rewritten in terms of the compactified ir-coordinate, the hypersurface 
equations for Q and W take the form 

x(l — x)d x F + 2F = RHS (2.11) 

where the right hand side is regular at . In order to deal with the degeneracy of this equation at x = 1, we rewrite 
(f2~TTj) in the form 

d(r 2 F) _ RHS 

d(r 2 ) " 2 [ ' 

and construct a centered finite difference approximation with respect to r 2 . Expressed in terms of the grid Xi — 
Xi-i + Ax, this leads to 

v _ ( Xj-\{1 - Xj) \ 2 y , Axjxj + x^i ~ 2x l x i _ 1 ) RHS 

''-Ua-si-ijJ ^(i-^-i) 2 2 [2A6) 

which enforces the correct asymptotic limit F\ x= i — RHS/2 when RHS is constant near J !+ . In practice, the 
variation of RHS implies that this limit is only enforced to first order accuracy when RHS is evaluated by the mid- 
point rule. This is consistent with global second order accuracy of Q and W when the numerical error is measured 
by an L2-norm over the hypersurface, but only first order accuracy can be expected for their values at ,y + . However, 
the asymptotic values of Q and W do not enter directly into the calculation of the waveform at y + . 

A. Waveforms at 

For technical simplicity, the theoretical derivation of the waveform at infinity is best presented in terms of an inverse 
surface-area coordinate £ — 1/r, where £ = at . In the resulting x^ — (u,£,x A ) conformal Bondi coordinates, 
the physical space-time metric g^ v has the conformal compactification g^ v = £ 2 g^ v , where g M „ is smooth at and, 
referring to (|2.ip . takes the form 24] 

g^dx^dx" = - (e 2p V£ 3 - h AB U A U B ) du 2 + 2e 2p dud£ - 2h AB U B dudx A + h AB dx A dx B . (2.14) 

As described in Q, the Bondi news function N(u,x A ) and the Newman- Penrose Weyl tensor component [43[ 

^1(u,x A ) = lim rip4 

r— >oo 

which describe the waveform are both determined by the asymptotic limit at of the tensor field 

% v = i(V M V„ - jg„ v V a V a )£. (2.15) 
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This limit is constructed from the leading coefficients in an expansion of the metric in powers of I. We thus write 

h AB = H AB + icAB + O(f). (2.16) 
Conditions on the asymptotic expansion of the remaining components of the metric follow from the Einstein equations: 

/3 = H + 0(l 2 ), (2.17) 

U A = L A + 2£e 2H H AB D B H + 0(£ 2 ) (2.18) 

and 

£ 2 V = D A L A + £(e 2H TZ/2 + D A D A e 2H ) + 0(£ 2 ), (2.19) 

where H and L are the asymptotic limits of /3 and U and where 1Z and Da are the 2-dimensional curvature scalar 
and covariant derivative associated with Hab- 

The expansion coefficients H, Hab, cab and L A (all functions of u and x A ) completely determine the radiation 
field. One can further specialize the Bondi coordinates to be inertial at J? + , i.e. have Minkowski form, in which 
case H — L A — 0, Hab — <1ab (the unit sphere metric) so that the radiation field is completely determined by cab- 
However, the characteristic extraction of the waveform is carried out in computational coordinates determined by the 
Cauchy data on the extraction worldtube so that this inertial simplification cannot be assumed. 

In order to compute the Bondi news function in the g^ v frame, it is necessary to determine the conformal factor uj 
relating Hab to a unit sphere metric Qab, i.e. to an inertial conformal Bondi frame [24[ satisfying 

Q AB =u 2 H AB . (2.20) 

(See [lij for a discussion of how the news in an arbitrary conformal frame is related to its expression in this inertial 
Bondi frame.) We can determine uj by solving the elliptic equation governing the conformal transformation of the 
curvature scalar (I2.7[) to a unit sphere geometry, 

K = 2(lu 2 + H AB D A D B loguj). (2.21) 

Equation (|2 . 2 1 1) need only be solved at the initial time. Then the geometrical properties of determines the time 
dependence of uj according to 

2h a d a logo; = -e- 2H D A L A , (2.22) 

where n a — g a/3 \7 p£ is the null vector tangent to the generators of J^ + . We use (|2.22l) to evolve ui along the generators 
of given a solution of (|2.21[) as initial condition. 

The news function N(u, x ) is first computed by the code in terms of the computational coordinates (u,x A ), as 
opposed to the inertial coordinates (u, y A ) on corresponding to an idealized distant observatory. The transfor- 
mation to inertial coordinates proceeds by introducing the conformally rescaled metric g^ u = uj 2 g^ v in which the 
cross-sections of J?^ have unit sphere geometry, in accord with (|2.20p . The rescaled null vector n" = u^ 1 n v is 
then the generator of the inertial time translation on i.e. h v d u — du- The inertial coordinates thus satisfy the 
propagation equations 

h v d v u = uj , n v d v y A = 0, (2.23) 

where h v d v — e~ 2H (d u + L a 8 x a) in terms of the computational coordinates. The inertial coordinates are obtained 
by integrating (I2.23[) . thus establishing a second pair of stereographic grid patches corresponding to y A . Then the 
news function is transformed into N(u 7 y A ). (More precisely, we should write N(u,y A ) = N(u,x A ) to distinguish the 
functional form of the news in the different coordinates but we forgo this complication of notation.) 

In addition, in order for the real and imaginary parts of N to correspond to the "plus" and "cross" polarization 
modes of a distant observatory, we need the proper choice of complex polarization vector Q' 3 , which in the inertial 
coordinates is related to the unit sphere metric on J'^ by Q AB — (Q A Q B + Q A Q B )/2. We fix the spin rotation 
freedom Q 13 -> e^Q? by requiring r^V^Q' 9 = 0(£l), so that the polarization frame is parallel propagated along the 
inertial time flow on .y + . This fixes the polarization modes determined by the real and imaginary parts of the news 
to correspond to those of inertial observers at ^ + . In order to carry this out in the computational frame we introduce 
the dyad decomposition H AB = (F A F B + F A F B )/2 where 




F A = l^r-^P 1 - tJ,l^-^- (2-24) 
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We then set — e iS oj 1 F 13 + Xh 13 , where F a := (0,0, F A ). The requirement of an inertial polarization frame, 
n v V ' v QP — 0(il), then determines the time dependence of the phase 8 according to 

2i(8 u + L A d A )S = D A L A + H AC F c ((d u + L B d B )F A - F B d B L A ). (2.25) 

The Bondi news now takes the explicit form 

N = ^e- 2i5 LJ- 2 e- 2H F A F B {(d u + £ L )c AB - ^c AB D c L c + 2ojD A [oj- 2 D B (uje 2H )}}, (2.26) 

where £ B denotes the Lie derivative with respect to L A . 

In the inertial Bondi coordinates, the expression for the news function reduces to the simple form 

N=^Q A Q B d u c AB . (2.27) 

However, the general form (I2.26[) must be used in the computational coordinates, which is challenging for maintaining 
accuracy because of the appearance of second angular derivatives of ui. 

Alternatively, the waveform can be obtained from the asymptotic value of the Weyl tensor. Asymptotic flatness 
implies that the Weyl tensor vanishes at J^ + , i.e. C^ up(y = 0(£) in the conformal Bondi frame (|2.14l) . This is 
the conformal space version of the peeling property of asymptotically flat spacetimes [23j]. Let (fi M , V 1 , m M ) be an 
orthonormal null tetrad such that n M = V^£ and l^d^ = di at . The radiation is described in this frame by the 
limit 



# :=-\\im\h»m v hPm a C fXUpa , (2.28) 
2 1^0 £ 



which in Newman-Penrose notation |43j corresponds to 

* = -(1/2)$. (2.29) 

The limit is independent of how the tetrad is extended off J> + . 
A major calculational result in [j| is that 

$=JnW(Vpt p -V„y| /+) (2.30) 



2 

where Y> a p is given by (|2.15p and where (|2.30[) is independent of the freedom 

rh u -> m v + Xh" . (2.31) 
In the same inertial polarization frame used in describing the news, 

* = ^- 3 e- 2iS h»F A F B (d^AB - d A t pB - f« B £ Aa + t'Xgt^ U+. (2.32) 

An explicit expression for ^ in terms of the asymptotic metric coefficients involves lengthy algebra which was carried 
out using a Maple script to write it in terms of 3 operators acting on the spin-weighted computational fields and to 
construct the final Fortran expression for 5*. 

In inertial Bondi coordinates, (|2.32j) reduces to the single term 

* = \Q A Q B d 2 u c AB = d 2 M^. (2.33) 

This is related to the expression for the news function in inertial Bondi coordinates by 

* = d u N. (2.34) 

However, as in the case of the news, the full expression for \P obtained from (|2 . 32[) must be used in the code. This 
introduces additional challenges to numerical accuracy due to the large number of terms and the appearance of third 
angular derivatives. 
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These difficulties can be appreciated by considering the linearized approximation, for which considerable simplifica- 
tion arises. To first order in a perturbation off a Minkowski background, the nonlinear expression (|2.32[) for ^ reduces 
to 

* = \dldiJ - \d u J - l^L - ig 2 (9Z + BL) + d u B 2 H. (2.35) 

Z Z Z a 

In the same approximation, the news function is given by 

N=^d u d e J+^d 2 (u + 2H). (2.36) 

The linearized Einstein equations imply that (|2.34p . i.e. = d u N , still holds in the linearized approximation. (In the 
nonlinear case, the derivative along the generators of is h v d v = e- 2H {d u + L A d A ) and (l233j) must be modified 
accordingly.) 

The linearized expressions (I2.35|) and (|2.36[) provide a starting point to compare the advantages between computing 
the radiation via the Weyl component "J or the news function N. The troublesome terms involve L, H and uj, which 
all vanish in inertial Bondi coordinates. One main difference is that \E' contains third order angular derivatives, e.g. 
the term 3 3 L, as opposed to second angular derivatives in the case of N. This means that smoothness of the numerical 
error is more crucial in the \f r approach. Balancing this, another main difference is that N contains the term S 2 w, 
which is a potential source of numerical error since U) must be propagated across the stereographic patch boundaries 
via (|2.22[) . Test comparisons of waveforms obtained via N and ^> are given in the next section. 



III. TESTS OF MODIFICATIONS TO THE STEREOGRAPHIC GRID 



The characteristic evolution carried out by the PITT code integrates the Bondi-Sachs equations by means of a finite 
difference approximation [2^, [26|. Stereographic coordinates x A — (q,p) are used to label the angles on the outgoing 
null cones. In the original code, two square stereographic patches were used, one centered about the North pole 
and the other about the South pole. In the new stereographic scheme introduced in [5[, the patches were modified 
to have circular boundaries located just past the equator, and angular dissipation was introduced to suppress the 
short wavelength noise introduced by interpatch interpolation. In addition, in the original code 9-derivatives were 
approximated by second order accurate finite difference approximations. In the present version used in this paper, 
the 3-derivatives have been increased to fourth order accuracy. Although the overall second order convergence rate 
of the PITT code remains unchanged, these changes are expected to lead to more accurate waveforms. 

There has been extensive testing of the accuracy of past versions of the code in [f| @, [H, [3l[ . Here we repeat some 
of the linear wave tests presented in [5J in order to demonstrate the improvement obtained by fourth order accurate 
angular derivatives. First, in order to verify that the new treatment of stereographic patches is capable of producing 
a fourth order accurate evolution, we carry out a test of wave propagation on the sphere based upon solutions to the 
2D wave equation 

-<9 2 $ + S9$ = 0, (3.1) 



where <i> = cos(ojt)Yi m , u> = yT{l + 1) and Yi m are spherical harmonics. For the case I = m = 2 we measure the 
convergence rate of the error. The simulations are run with n + 1 grid points along the axes of each patch, with the 
grid sizes ranging from n = 80 to n = 240. For a given grid size, we use the norm to measure the error 

£(*^) — 1 1 ^numeric ^ analy tic \ \ oo (^-^) 

for the circular patches in each hemisphere. We measure the convergence rate for £ ($) at a given time t, for two 
consecutive grid sizes ri\ and ri2, by 

_ log 2 (£($) n2 /£($) ni ) 
log 2 (ni/n 2 ) 

Convergence rates for the derivatives are measured analogously. 

Excellent 4th order convergence of £($) was obtained. It is more important and challenging for assessing waveform 
extraction error to measure the error in the derivatives 9$, <3 2 $, and <3 3 <I>, since second angular derivatives enter in 
the computation of the Bondi news and third angular derivatives enter into the computation of 5'. The convergence 
rates, measured with the L m norm over the North patch, are shown in Table U based upon the grid sizes (ni,ri2) = 
(80, 120), (120, 160), (160, 200), (200, 240). 
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TABLE I: Convergence rates for errors in 0$, 3 2 $ and 3 3 $ 



S/n 


m = 80 m = 120 m = 160 m = 200 


£(a 2 $) 

£(3 3 $) 


4.04 4.11 4.35 4.85 
4.07 4.24 4.80 3.95 
3.98 3.95 3.92 3.86 



For the coarser grids, good 4th order convergence is apparent for all the derivatives. As the grids are refined, the 
error eventually approaches (double precision) roundoff error and convergence becomes a moot question. 

Next we compare the accuracy of waveform extraction by computing the news function N and the Weyl tensor 
component \& in the test problem considered in which is based upon a periodic, linearized gravitational wave on a 
Minkowski background (see Sec. 4.3 of [45j]). The linearized wave is expressed in Bondi-Sachs coordinates so that it 
allows direct measurement of the numerical error. The wave has period T = tt and (I = 2, m = 0) spherical harmonic 
dependence, with the maximum value of J rs 1CP 6 . The data provided by the linearized solution at the extraction 
worldtube was propagated to by the characteristic code, where the waveform was computed and compared to its 
analytic value. The computational error in the waveform was measured with the L2 norm over the North patch using 
the n = 100 grid. 

Figure [T] compares plots of the error in the real part of the news function TV (computed on the North patch) for 
the 2nd and 4th order accurate angular derivatives. The plots show roughly one order of magnitude improvement in 
accuracy for the n — 100 grid. The corresponding plots of the error in the waveform measured by the Weyl component 
'J show again roughly one order of magnitude improvement in accuracy. Further improvement in accuracy might be 
obtained by also increasing the radial derivatives to fourth order approximations but this could entail nonlinear 
complications which could affect the numerical stability of the evolution algorithm [42| . 




FIG. 1: Plots of the L2 errors £(N) vs t in the real part of the news function extracted in a linearized gravitational wave test. 
The plots compare the errors obtained using the 2 nd and 4 th order accurate angular derivatives on an n = 100 grid. The fourth 
order method reduces the error by an order of magnitude. The time variation of the error matches the period of the wave. 



In accord with ()2.34|) . the computation of the Weyl component ^ yields an alternative numerical value for the news 

ru 

Nq, = N\ U=Q + / *du, (3.4) 
Jo 

where N = in the analytic problem. Figure [2] compares these two extraction methods in terms of the errors in N 
and for the linearized wave test when using 4th order accurate angular derivatives. The plots show that the two 
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methods are competitive although the error in Nq, is slightly smaller in this case. 




FIG. 2: Comparison plots of the L2 errors £(N) and £(N\j,) vs t in the news function computed directly and via the Weyl 
tensor for the linearized gravitational wave test. The results were obtained using 4 th order accurate angular derivatives on an 
n — 100 grid. The two methods are competitive although £(N<&) is slightly smaller in this case. 



IV. COMPUTATIONAL INTERFACE 



We have designed an interface that takes Cartesian grid data from a Cauchy evolution and converts it into boundary 
data for a characteristic evolution on a spherical grid extending to J" + . We treat each component g^it^x 1 ) of the 
Cauchy metric as a scalar function in the x % Cartesian coordinates which are used in the 3 + 1 evolution. In order to 
make the interface as flexible as possible for use as a community tool for waveform extraction, we have based it upon 
a spectral decomposition of the Cauchy data in the region between two world tubes or radii R = R\ and R = R2, 
where R = yj8ijX l x^ is the Cartesian coordinate radius. Then at a given time t, we decompose g^t.x 1 ) in terms 
of Chebyshev polynomials of the second kind Uk{R) and spherical harmonics Yi m (6,<p), where (9,4>) are related to 
x l I R in the standard way. The Chebyshev polynomials are conventionally defined as functions C/fc(r) on the interval 
— 1 < r < 1. Here we map them to the interval R\ < R < Ri by the transformation 



Ri — R\ 



Thus, for Ri < R < R 2 , we expand 



g^{t,x l ) = y2 C Mklrn]{t)U k (R)Y lm (6,<t>). (4.1) 



klm 



For the application to waveform extraction given in this paper, we choose I < IMax, where IMax — 6, and k < kMax, 
where kMax — 6. These values should be considered tentative and further experimentation is warranted to optimize 
accuracy. In tests of binary black holes with mass M we use a relatively small range R2 — R\ = 10 M and a larger 
value of kMax would certainly be needed if the range were expanded. Also, while IMax = 6 might be sufficient for 
extraction at Re = 100M, a larger value might give improved results at Re = 20M. 

The coefficients C^kim] allow us to reconstruct a spherical harmonic decomposition of each component of the 
Cauchy metric on the extraction worldtube R — Re, i.e. 



u[klm] 

(t)U k (R E ). 

k 



(4.2) 
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This decomposition is carried out at a sequence of Cauchy time steps t n = to + nAt, where At is chosen to be much 
smaller than the physical time scales in the problem but, for purpose of economy, larger than the time step used for 
the Cauchy evolution. At each time step, the spectral coefficients are determined by a least squares fit to the Cauchy 
metric. 

The extraction module also requires the derivatives dtg^v and dng^u at the extraction worldtube. The i?-dcrivative 
is obtained analytically, at each time level t n , by differentiation of the Chebyshev polynomials. In one option, the 
finite difference option, the t-derivative is constructed by a fourth-order accurate finite difference approximation 
based upon the sequence of Cauchy times t = t n . In a second option, the fast-Fourier-transform option, we modify 
the Cauchy data by filtering each mode f n = Cki m {t n ) to remove high-frequency noise (both numerical noise and 
high-frequency gauge waves). The filter works as follows. Let /„ be the original data (n = 0, ...,N — 1), and 
9n — fn ~ a(nAt) — b(n At) 2 . The coefficients a and b are fixed by requiring that g^ = go (where g^ is extrapolated 
from gN-i and gN-2) and go — gjv-i = .92 — 9i, i-e. the one sided derivatives taken at n — agree. This guarantees 
continuity of g and its first derivative when periodically extended. We then perform a fast Fourier transform on gi, 
truncate the transform at high frequencies, and perform an inverse Fourier transform to obtain a filtered G n and 
optionally, the inverse transform of iu) gi to obtain a smooth time derivative of G n . We then construct the filtered 
mode Ckim(t n ) = G n + a(nAt) + b(nAt) 2 , as well as its time derivative. 

The stereographic coordinates x A — (q,p) used to label the outgoing null rays in the Bondi metric are matched 
to the spherical coordinates {9, <f) induced by the Cartesian Cauchy coordinates on the extraction worldtube by a 
standard transformation, using the conventions in |40j |. The value of the surface-area coordinate r in the Bondi- 
Sachs metric is obtained on the extraction worldtube from the 2-determinant of the Cartesian metric on the surfaces 
t — t n ,R — Re- As a result r E (t n ,q,p) := r\n = R E 7^ const on the extraction worldtube. In order to deal with this 
complication, the transformation from Cartesian coordinates (t,x l ) to Bondi-Sachs coordinates (u, r, x A ) is carried 
out via an intermediate Sachs coordinate system (u,X,x a ) [461 ] where A is an affine parameter along the outgoing 
null rays. The affine freedom allows us to set A = on the extraction worldtube. After carrying out the Jacobian 
transformation from (t, x l ) to (u,X 7 x A ), the Cartesian metric and its first derivatives at the extraction worldtube 
provide a first order Taylor expansion in A (about A = 0) of the null metric in Sachs coordinates. The corresponding 
Taylor expansion of the metric in Bondi-Sachs coordinates then follows from the computed values of te and d\r at 
A = 0, which are obtained from the 2-determinant of the Cartesian metric Q. 

This allows us to build a grid based upon the characteristic coordinates [x, q,p), with compactified radial coordinate 
x given by (|2.9p . The grid values Xi = a^-i + Ax, 1 < i < n Xl are adjusted so that x\ < xe '■= x\r = r e and x ni = 1. 
The characteristic time levels u n = + Au are chosen to coincide with the Cauchy times t n on the extraction 

worldtube by choosing (u — t)\n = R E = 0. 

In the previous version of the extraction module, the first order Taylor expansion for the Bondi metric was used 
to fill the gridpoints neighboring the extraction worldtube and thus initiate the radial integration of the hypersurface 
equations (|2.3|) - (|2.6[) . However, the hypersurface equations require only 6 (real) integration constants, which can be 
supplied by their values at R = Re- Using the Taylor expansion to fill the neighboring gridpoints leads to a potential 
inconsistency between the Bondi metric supplied by the Cauchy evolution and the radial derivatives determined 
by the characteristic hypersurface and evolution equations. In particular, we have found that such inconsistencies 
arising from error in the Cauchy data degrade the convergence rate of the characteristic extraction module. Because 
convergence of the extraction module is an important test of its reliability, we proceed here in a different manner 
which decouples the Cauchy and characteristic extraction errors. 

In the previous version of the extraction module, the Taylor expansions were also applied to the auxiliary variables 
v = 8 J, B = 5/3 and k = dK by applying the 3-operator to the Taylor expansions of the main variables. This was 
a complicated process because the S operator intrinsic to the A = extraction worldtube is not the same as the 3 
operator intrinsic to the r = const Bondi spheres (as they differ by radial derivatives). In the process, several bugs 
were introduced in the radial start-up scheme. The present version of the extraction module streamlines the start-up 
of the auxiliary variables by avoiding the use of Taylor expansions. 

In this new approach, the hypersurface equations are integrated purely in terms of the values /3e, Qe, Ue and We 
of the hypersurface variables on the extraction worldtube which are supplied by the Cauchy data. A mask is set up 
to identify those radial grid points i < B (referred to as U B points") for which Xi — xe < Ax. These grid points are 
"passive" points which do not directly enter in the evolution. Values of the hypersurface variables are assigned at 
the first active points i = B + 1 (referred to as ll B + 1 points") in the following manner, assuming that the values 
Je and Jb+i of the evolution variable J are known, as well as the values vb+i and ks+i of the auxiliary variables. 
(We address the latter assumption below in describing the start-up of the evolution algorithm.) Proceeding in the 
hierarchical order of the hypersurface equations, we first use (|2 . 3() to determine /3b+i according to 

p B+ i=PE + Nf)[J](rB+i-r E ). (4.3) 

Because Np[J] only involves J and d r J it may be evaluated at the mid-point between xe and xb+i so that the 
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resulting error in (3b+i is 0(Ax 3 ). This also determines the auxiliary variable B — 9/3 at the B + 1 points provided 
the B + 1 points on the neighboring rays have the same grid value However, in the case of an irregularly shaped 
extraction worldtube, there can be exceptions where this neighboring ray is a B point. As a result, in cases where 
the B points lie close to the boundary of the masked region they can couple to the B+l points on neighboring rays 
through the 3 operator. For this reason, we also update B points by the same scheme used for the B + l points. (If a 
B point is within a small tolerance of the world-tube, we instead just copy the world-tube value rather than risk an 
ill-conditioned algorithm.) In this way, the start-up value of the auxiliary variable Bb+i is determined in all cases. 

Next in the hierarchy of hypersurface equations, we determine Qb+i in similar fashion. However, Nq[J, /?] involves 
B = 9/3 which cannot be determined on the extraction worldtube from the values of Pe (because of the angular 
variation of te discussed above). Consequently, in order to start up the Q-intcgration we evaluate Nq at xb+i, where 
Bb+i is known. This results in an 0(Ax 2 ) error in the value of Qb+i- Similar considerations apply to the start-up 
of the U and W integrations. As a result, the start-up leads to an overall 0(Ax 2 ) error in values at xb+i, which 
is consistent with the global 0(Ax 2 ) error resulting from the remaining integration from xb+i to This radial 

march to y + proceeds in the same way as described in (2(| EH, EH to determine all variables on the hypersurface. 

Having completed the radial march on the hypersurface at time un-i, the start-up of the integration scheme on 
u n — u„_i + Au begins with the determination of JB+i{u n ) from the worldtube data Je, /3e, Qe, Ue and We on 



un and the fields already determined on un-i- We determine Js+i{u n ) using a null parallelogram algorithm 47 1. 
The evolution equation (|2.8p for J can be rewritten as 

2d u d r <I> - d r (Ad r <S>) = RHS (4.4) 

where $ = rJ and A = V/r = 1 + rW. This can be integrated over the null parallelogram in the (u,r) subspace 
bounded by the u n and it n _i hypersurfaces and by two ingoing characteristics. For constant A, the ingoing charac- 
teristics satisfy r — (Au/2) = const. As depicted in Fig. [3j by choosing one ingoing characteristic to pass through rg 
on the u n hypersurface and the other to pass through Tb+i on the midpoint between the u n and hypersurfaces, 
we obtain the integral approximation 

^/ ^ \ i \ RHS(r+ — r )Au 

*(«n,r_) = $(u n ,r E ) +*(u„_i, r+) - $(u„_i,r ) + ^- - . (4.5) 

Here the corners of the null parallelogram are located at r E , T± = tb+i ± (AAu/A) and tq = r E + (AAu/2); and the 
center of the null parallelogram is located at r c = (l/2)(r^ + r+). This determines the start-up value $(tt n ,rB+i) 
through the second order accurate interpolation 

3 ( j = [*(un,rB + i)-*(un,r g )]r._ 

Tb+i - r E 

Using the worldtube data and field values on m„_i, all other quantities can be approximated consistent with second 
order accuracy except for a term in RHS which is proportional to d u J. This term is treated to the required accuracy 
by a two-step Crank-Nicholson iteration, as is done in the main evolution scheme described in [26| . This leads to 
a value of $(«„, Tb+i), and thus J(u n , f E +i), with 0{AxAu 2 ) error. As in the case of the hypersurface equation, 
we also use this algorithm to update J at the B points to assure that the auxiliary variables vb+\ and ks+i can be 
determined by application of the 3 operator. Now the radial march continues to the B + 2 points by a similar process. 



A. Convergence measurements 



In tests of the waveform and other variables obtained from a binary black hole evolution there are no exact values 
available for measuring error so that convergence rates cannot be obtained by use of (|3.3I) . Instead, we obtain 
Cauchy convergence rates by using measurements obtained with three different gridsizes. For grids in the ratio 
A3 = x^2 = X 2 A\ the Cauchy convergence rate of a measured quantity F is given 

n = l0g 2 ((^3-F 2 )/(F 2 -Fl)) 

log 2 X 

For quantities that approach the continuum value Fq as F = Fo + GA 2 , (14.71) gives a convergence rate of 2 when G is a 
smooth function independent of gridsize. In the main part of the characteristic evolution algorithm, G is determined 
by the second derivatives of the evolution variables. However, a stochastic grid-dependent source of second order error 
occurs in the start-up algorithm due to the location of the B + l points. The separation xb+i -i£ of this point from 
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FIG. 3: The start-up diagram in the (it, r) subspace for the B + 1 points (shaded squares). On the left, the extraction worldtube 
with fixed Cartesian radius Re moves with respect to the null grid. The null parallelogram for the start-up algorithm is bounded 
by the two outgoing characteristics at retarded times w„_i and u n and the two ingoing characteristics indicated by dashed 
lines. The labels for the radial null coordinate r are indicated at the four corners (shaded circles). 



the extraction worldtube can vary discontinuously under a small change in gridsize, i.e. xb+i -ib = (l + e)Aa;, where 
e is a random number, < e < 1. This random separation enters into the second order accurate approximations made 
in the start-up algorithm. The approach to the continuum value has the form F = F + (G + Ge)A 2 . Consequently, 
the stochastic part of the second order error can obscure the convergence rate determined by (14.71) if G is comparable 
in size to G. The only way to ensure clean convergence rates would be to implement a third order accurate start-up 
algorithm, which would involve a considerable amount of work. Fortunately, this source of error does not appear to be 
significant in the tests we have carried out. All the main variables exhibit second order convergence when measured 
at a finite radius for the results of the binary black hole inspiral presented in Sec. [V] However, some asymptotic 
quantities at ^ + display convergence rates intermediate between first and second order, for reasons discussed further 
in Sec. ED 

B. Constraints on the time step 

Domain of dependence considerations place a constraint between the characteristic time step Am and the size of 
the characteristic grid analogous to the CFL condition for the Cauchy evolution. For a rough estimate, consider the 
Minkowski space case with the conformally rescaled metric 



(1 - x) 2 2 

ds 2 = ~ p2 du 2 - —dudx + qAsdx A dx B (4.8) 
R E Re 



in the compactified stereographic coordinates (u, x, q,p) used in the code, for which the unit sphere metric takes the 
form 



q AB dx A dx B = —4- 2 (dp 2 + dq 2 ). (4.9) 
1 + V + q z 



The past light cone is determined by 



du —dx — y/ 'dx 2 + (1 — x) 2 qAgdx A dx B 

Re~ = (T^P ' 
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The restriction on the characteristic time step arising from domain of dependence considerations is strongest at the 
inner boundary, where x = 1/2 (since te = Re in the Minkowski case); and it is also strongest at the equator, where 
p 2 + q 2 = 1. At these points 

\du\ 



dx + yjdx 2 + (1/4) (dp 2 + dq 2 ). (4.11) 



4i? s 

For typical characteristic grid parameters, Ap = Aq = Ax/4, the resulting restriction is 

\Au' 



Re 



< 8KAx (4.12) 



where K ~ 1 depends upon the details of the finite difference stencil. This restriction is strongest for a small extraction 
radius. The characteristic code monitors the corresponding restriction on Am determined by the curved space version 
of the compactified Bond-Sachs metric. For a Cauchy simulation of binary black holes with mass M with timestep 
At = M/32 (sufficient to describe the frequencies typical of a binary system), (I4.12[) leads to 

M 

< KAx, (4.13) 



256i? B 



for the choice of characteristic timestep Au = At. The corresponding number of radial gridpoints must roughly satisfy 
n x < 128.Re/M. This places no limit of practical concern on the resolution of the characteristic evolution even for 
the small extraction radius Re = 20M. Thus, for purposes of CCE, there are no demanding CFL restrictions on the 
characteristic grid and timestep. 



C. Initial characteristic data 



The initial data for characteristic evolution consist of the values of J on the hypersurface u = Tq. One way of 
attempting to suppress incoming radiation in this data is to set the Newman-Penrose Weyl tensor component ipo = 
on the initial null hypersurface. For a perturbation of the Schwarzschild metric, this condition implies no incoming 
radiation in the linearized approximation. The condition that ?/>o vanish involves two-radial derivatives of J, which 
in the compactified coordinate I = 1/r takes the simple linearized form d 2 J —0. Translated into the computational 
coordinate x = 1/(1 + Re^), we choose the solution 

J= J\ x (l-x)x E 
(1 - x E )x 

which provides continuity of J with its value determined by Cauchy data at the extraction worldtube. Since this 
choice of J also vanishes at infinity, the initial slice of J" + has unit sphere geometry and equation (12.21)) for the 
conformal factor has the simple solution co = 1. 



V. BINARY BLACK HOLE MEASUREMENTS AND WAVEFORMS 



Here we present test results of waveform extraction from the inspiral and merger of two equal-mass, non-spinning 
black holes. For the Cauchy evolution we used the LazEv code (la , IHij along with the Cactus framework [5(| and 
Carpet [5l[ mesh refinement driver. LazEV is a 8th order accurate finite difference code based upon the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formulation [l9l [20j of Einstein's equations, which deals with the internal singu- 
larities by the moving puncture approach (48l |52|. Our simulation used 9 levels of refinement with finest resolution 
of h = Af/80.64, and outer Cauchy boundary at 400Af . The initial data consisted of a close quasicircular black-hole 
binary with orbital frequency MVt = 0.050, leading to more than a complete orbit before merger (See (53|). We output 
the metric data on the extraction worldtube every At = M/20. 

In the Cauchy evolution, we extract ip4 on spheres of Cartesian radius R/M = 50, 60, • • • , 100 and decompose in 
spin-weighted (/, m) spherical harmonic modes. We use a perturbative formula 55] to extrapolate the perturbative 
waveform Rip 4 to R — oo, 

lim [Ri> l 4 m (R,t)] = r^{ m (r,t) - ( * ~ ^ + 2) / dt^\r,r)dr + 0(r 2 ), (5.1) 

R^oo 2 Jr. 
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where r is the areal radius corresponding to the Cauchy extraction radius R. The extrapolation of the perturbative 
waveform to infinity removes cumulative phase error which otherwise would be introduced by redshift effects. 

We present results for the characteristic extracted waveform either in terms of related to the Bondi news by 
\& = d u N, or, when comparing to the perturbative waveform, in terms of the Newman-Penrose component "J 4 = — 2\E'. 
For illustrative purposes, we concentrate on the dominant (I = 2,m = 2) and sub-dominant (I = 4, m = 4) spherical 
harmonic modes. 

The highest resolution black hole waveform extraction test was run with the following characteristic grid specifi- 
cations: angular gridpoints n q = n p = 200, radial gridpoints n x = 224. For convergence tests, we also used grids 
n q — n p — 100, n x — 112 and n q = n p = 50, n x = 56, so that the grid sizes were in the ratio \ = 2. We re- 
fer to these as the n = 200, 100 and 50 grids, respectively, The characteristic time steps used for these grids were 
At = M/20 (n = 200), At = AT/10 (n = 100) and At = AT/5 (n = 50). The characteristic extraction was carried out 
using worldtube radii R E = 20Af, 50M and lOOAf . 

The Pitt null code was run on stereographic patches with circular boundaries using the auxiliary variables (|2.10p to 
eliminate any second derivatives in the angular directions and using 4th order accurate angular derivatives. Angular 
dissipation was added with the coefficients e x = 10~ 3 , e u = 10~ 4 , eq = e\y = 10~ 6 , in the notation of 

The best accuracy was obtained using the fast Fourier transform (FFT) option to obtain time derivatives of the 
worldtube data, as described in Sec. IIV1 For strong signals, e.g. the dominant (I = 2, m = 2) spherical harmonic 
mode, the finite difference (FD) and FFT options are in good agreement. However, for weak signals, e.g. the 
(/ = 4, m = 4) mode, the FD option can generate noticeable high frequency error. See Fig. 2] for a comparison of 
waveforms computed with the two options. One likely source of the high frequency error with the FD option is the 
stochastic error introduced at each time level t n on the extraction worldtube by the least squares fit of the Cauchy 
data to the spectral expansion. In the FD option, this error is amplified when taking the time derivatives necessary 
to compute the waveform and becomes more prominent for short characteristic timesteps. It also becomes more 
prominent as the extraction radius is increased, in which case the extracted worldtube data is smaller. Similar high 
frequency noise is apparent in the worldtube data so that this error cannot be removed by refining the characteristic 
grid. However, the filtering intrinsic to the FFT option is effective in reducing this error. The remaining results 
reported in this paper were obtained with the FFT option. 
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FIG. 4: Comparison of the FD and FFT options for the (2, 2) (up) and (4, 4) (down) spherical harmonic modes of the real 
part of the characteristic waveform $ obtained with the n — 200 grid and extraction radius Re = 20M. The two options give 
comparable results for the (2,2) mode but for the (4,4), which is an order of magnitude smaller, high frequency error in the 
world tube data is noticeable in the waveform extracted with the FD option. 

Convergence rates were measured for the (I = 2, m = 2) spherical harmonic mode, which is the dominant mode in 
the waveform. Table |H] gives the rates for the evolution variables measured on a sphere at Bondi radius r = 80M 
obtained with a small extraction radius Re = 2QM at a time corresponding to the peak of the signal (t « 200A/). 
The rates are given for the real and imaginary part of the variables. All quantities are very close to second order 
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convergent, including J )X , which is the term which determines the waveform after transforming to inertial Bondi 
coordinates according to (|2.33[) . 

TABLE II: Convergence rates of the (I = 2, m = 2) spherical harmonic mode on the sphere r — 80M for the metric variables 
measured at retarded time u ~ 200M near the peak of the signal. The rates are given for the real and imaginary part of the 
variables. The extraction radius was R — 20M. The results show that the evolution variables all display clean second order 
convergence. 



Variable 


RateR e 


Rateim 


P 


2.01 


2.01 


J 


2.23 


2.01 


J ,X 


2.03 


2.33 


Q 


2.02 


2.04 


u 


1.99 


1.96 


w 


1.97 


2.00 



Table UTTl gives the corresponding convergence rates for these evolution variables measured at again at the time 
corresponding to the peak of the signal and with extraction radius Re = 20M. In this case Q and W show deviation 
from second order convergence, consistent with the asymptotic error analysis presented in Sec. Hi] in relation to (|2.13|) . 
We also see that the derivative J jX deviates from second order convergence, which indicates a need for more accurate 
finite difference approximations near There are several places in the present code where one-sided difference 

approximations are used for derivatives at ^ + . These convergence rates at the peak of the signal are representative 
of the rates over the entire run. This is illustrated in Fig. [5] which plots the rescaled errors of Re J and ImQ versus 
time at 

TABLE III: Convergence rates of the (I = 2, m = 2) mode for the metric variables measured near the peak of the signal at 
J^ + , with an extraction radius R = 20M. As expected from the analysis in Sec. [TT] some asymptotic quantities only display 
first order convergence. 



Variable 


RateRe 


Rateim 


P 


2.01 


2.01 


J 


1.80 


2.18 


J ,X 


1.23 


1.20 


Q 


1.33 


1.19 


u 


1.99 


1.96 


w 


1.55 


1.50 



Table IIVI gives the corresponding convergence rates for the waveform as measured by the Bondi news N and the 
Weyl component again at a time corresponding to the peak of the signal and with extraction radius Re = 20M. 
We also show the convergence rate of the inertial time derivative d u N calculated directly from finite differencing N. 
All show roughly first order convergence. The rate for d u N is slightly better than that for \&, although \& = d u N in 
the continuum limit. The convergence rates of these quantities are affected by two chief factors: (i) the large number 
of terms involved in their calculation and (ii) their dependence on radial derivatives of the evolution quantities at J? + . 
In all cases, one-sided approximations are used in several places to compute these radial derivatives. This is already 
apparent in the convergence rate for J jX shown in Table IIIII 

Surface plots of the Bondi news N and Weyl component VP, measured near the peak of the signal with an extraction 
radius Re — 20M, are shown in Fig's [5] and [7J Both figures display smooth angular dependence, showing that the 
angular dissipation is effective at removing short wavelength angular noise. In particular, there are no "spikes" near 
the equatorial patch boundary arising from interpatch interpolation. The main error in the waveform originates from 
intrinsically time dependent error in the data on the extraction worldtubc. 

The time dependence of the real part of the characteristic extracted waveform and its comparison to the pertur- 
bative waveform are shown in Fig's [8] and |9] Figure [8] shows excellent agreement between these waveforms for the 
dominant (I = 2, m — 2) mode, when both are extracted at R = 50M. The insets show that this agreement exists 
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FIG. 5: Convergence plots of the asymptotic limits at of the (2, 2) spherical harmonic modes of ReJ and ImQ obtained 
with resolutions n = 50, n = 100 and n = 200 with an extraction radius Re = 20M. The plots for ReJ are rescaled for 
2nd order convergence (upper plot), while the plots for ImQ (lower plot) are rescaled for 1st order convergence. The rescaled 
differences show that convergence rates at the peak of the signal given in Table 11111 are representative of the rates over the entire 
run. 



TABLE IV: Convergence rates of the (2, 2) spherical harmonic mode for the Bondi news N, d u N obtained by finite difference, 
and the Weyl component all measured near the peak of the signal with an extraction radius Re = 20M. 



Variable 


RateRe 


Ratei m 


N 


1.59 


1.56 


d u N 


1.57 


1.55 




1.16 


1.14 



in the early stages, when the amplitude is small, and persists throughout the final ringdown. Note that the pertur- 
bative extrapolation formula (|5.ip is essential to obtain this excellent phase agreement between the perturbative and 
characteristic waveforms. 

Figure |H] compares the characteristic and perturbative waveforms for the (I = 4, m = 4) mode. In this case, the 
perturbative waveform is again extracted at R — 50M but the characteristic waveform is extracted at 20M to reduce 
the high frequency noise discussed previously. This high frequency noise can also be reduced by choosing a larger 
timestep for the characteristic evolution, again indicating that it arises from time derivatives of the stochastic error 
introduced in the worldtube data by the least squares fit. The characteristic and perturbative waveforms again show 
excellent agreement At very early times t/M « 15, the characteristic waveform shows an anomalous feature which 
can be attributed to "junk" radiation content in the initial data in the vicinity of Re = 20 A/. In addition, as shown 
in the insert, there is another anomalous feature, which is not understood, in the time interval about t/M = 90. This 
feature is also evident in the extracted Cauchy data at Re = 20M. Possible sources for this feature are gauge modes 
excited in the interior region or numerical effects from the adaptive mesh refinement used in the Cauchy evolution. 
A better explanation would require further runs. 



VI. RICHARDSON EXTRAPOLATION AND CONVERGENCE OF THE WAVEFORM 



The clean first order convergence results for the news N and Weyl component \& allows us to apply Richardson 
extrapolation to obtain higher order accuracy waveforms. We apply the results from the three different resolutions 
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FIG. 6: Surface plot in the North (q,p) stereographic patch of the real part of the Bondi news N measured at the peak of the 
wave with an extraction radius Re = 20M. The equatorial patch boundary corresponds to the circle p 2 + q 2 = 1. The smooth 
angular dependence near the equator shows that angular dissipation is effective at removing short wavelength noise arising 
from interpatch interpolation. 



n = (50, 100, 200), with grid spacing (4A, 2A, A) respectively, to obtain a third order accurate waveform as follows. 
The truncation error in a quantity F can be represented by a power series 

F(A) = F + F'A + If" A 2 + 0(A 3 ). 

We write Fi = /(A), F 2 = F(2A) F 4 = F(4A). Then the extrapolated value 

8 1 
Fp: = —F-] — 2_ro H — Fa 
3 3 

is 3rd order accurate, i.e. 

F E = F + +0(A 3 ). 

In practice this can be confirmed as follows. Let Fj = 2F\ — F2 and F/j = 2F 2 — F4 be the second order accurate 
waveforms obtained using data from two resolutions. Then F/j — Fe = 4(F/ — Fe) + 0(A 3 ), i.e. 

\{F U - F E ) = F/ - F E (6.1) 

if we neglect the 0(A 3 ) error, i.e. if we approximate the exact value Fq by the third order accurate approximation 
F E . 

Using the corresponding notation (Ng, Nj, Njj) for the news and (^>e, ^1,^11) for the Weyl component, we can 
check the validity of applying Richardson extrapolation to the waveform. Figure [TU1 and Figure [TT] graph the rescaled 
errors of the real and imaginary parts of Nj(t) — iV^(t) and k(Nu(t) — Nsit)) and Figure [12] graphs the corresponding 
rescaled errors in ^(t). In both cases, (I6.1[) is confirmed. 

These results validate the use of Richardson extrapolation to obtain third order accurate waveforms Ne and e- 
We can use Ne and as fiducial exact values and estimate the truncation error in the numerical waveforms by 
comparing them with the second order accurate approximates A/j and ^j. Thus the truncation errors in the news 
iVand Weyl component "3/ are conservatively given by 



5N = Ni — N E — 0(A 2 ) 



(6.2) 
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FIG. 7: Surface plot in the North (q,p) stereographic patch of the real part of the Weyl component ^ measured at the peak 
of the wave with an extraction radius Re = 20M. The smooth angular dependence near the equator p 2 + q 2 — 1 shows that 
angular dissipation is effective at removing short wavelength noise arising from interpatch interpolation. 




FIG. 8: Comparison of the (2,2) dominant spherical harmonic mode for ^4 (characteristic) and ripi (perturbative Cauchy), 
both extracted at R = 50M. The insets show that the excellent agreement extends to the early stages and the final ringdown 
when the amplitude is small. 
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FIG. 9: Comparison of the (4, 4) sub-dominant mode for ^4 (characteristic) extracted at Re = 20M and ripi (perturbative 
Cauchy) extracted at R = 50M. There is good agreement in the strong amplitude regime of the wave t/M > 120. At early 
times t/M ~ 15, the characteristic waveform exhibits effects of "junk" radiation in the initial data near Re = 20M. In addition, 
the insert magnifies an anomalous feature, which is not fully understood, in the interval about t/M = 90. 




FIG. 10: Plots confirming the validity of Richardson extrapolation to obtain higher order accuracy for the real and imaginary 
parts of the dominant (2, 2) spherical harmonic mode of the news N(t). The rescaled errors show that Ni and Nn are second 
order accurate in accord with (|6.ip . 
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FIG. 11: Plots confirming the validity of Richardson extrapolation to obtain higher order accuracy for the sub-dominant (4, 4) 
spherical harmonic mode of the news N(t). The rescaled errors show again that Ni and Nn are second order accurate in 
accord with (16.11). 



and 



6V = = 0(A 2 ). 



(6.3) 



Figure [13] plots the differences between the dominant (I = 2, m = 2) mode of the Richardson extrapolated waveform 
NE(t) obtained with extraction radii Re = 20M, Re = 50M and Re = 100M. In the plot, the Re = 20M waveform 
begins at t = and the other waveforms have been shifted backwards in time so that all three are in phase at the 
peak of the wave. Two sources of extraneous "junk" radiation can be seen in the figure. One arises from a mismatch 
between the initial characteristic and Cauchy data. The initial characteristic data ^0=0 (see Sec. IIV C[) implies 
the absence of initial radiation content on the assumption that the geometry of the initial null hypersurface is close 
to Schwarzschild. This assumption becomes valid as the extraction radius becomes large and the exterior Cauchy 
data can be approximated by Schwarzschild data. Thus this mismatch is largest for extraction at Re = 20M. This 
results in a noticeable difference at very early times between extraction at Re — 20M and the other two radii. After 
t/M « 100, the three waveforms are in good agreement with their relative differences less than 0.6% at the peak of 
the wave. 

The second source of "junk" radiation apparent in Figure [T3] arises from the choice of conformally flat initial Cauchy 
data. This arises for all three extraction radii and accounts for the double hump in the news function in the interval 
from t/M = to t/M = 50. 

It is of interest to measure the difference 



(6.4) 



between the extracted waveform using the perturbative extrapolation formula (15.11) and the Richardson extrapolated 
characteristic waveform, measured in accord with the normalization conventions indicated in (|2.29|) . Figure [14] plots 
the real part of the (2,2) spherical harmonic component of 8ipi{t) 1 compared with the corresponding component of 
Rety . The peak amplitude of Sip^t) is approximately 1% of the peak amplitude of 'F, which provides an estimate of 
the difference between perturbative and characteristic extraction. 

Figure [15] plots the phase difference S$ between the (2, 2) spherical harmonic components of r04(i)/2 and ^(t), i.e 



(5$ = - $[</> 4 ], 



(6.5) 



where e.g. <F = |\l/|e z4 '[*l. The phase difference is less than 0.05 radians in the interval 40M < t < 250M beginning 
after the initial burst of junk radiation and extending into the late ringdown. The phase errors become large at late 
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FIG. 12: Plots confirming the validity of Richardson extrapolation to obtain higher order accuracy for the real and imaginary 
parts of the (2, 2) spherical harmonic mode of the waveform >&(£). The rescaled errors show that and vp/j are second order 
accurate in accord with (|6.ip . Note that the second order error in ^ contains more high frequency noise than TV shown in 
Fig. [lOl 
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FIG. 13: Plots of the (2,2) mode of ReN obtained for extraction radii R E = 20M, 50M, and 100M. The R E = 50M and 
Re = 100M waveforms have been shifted backward in time so that they are in phase at the peak of the wave. The noticeable 
difference in the interval from t/M = to t/M = 100 between the Re = 20M waveform and the other two results from a 
mismatch between the initial characteristic and Cauchy data, which decreases with large extraction radii. For the waveforms 
extracted at all three radii, the double hump in the interval from t/M — to t/M = 100 results from non-trivial "junk" 
radiation in the initial Cauchy data. The three waveforms are in good agreement in the inspiral and merger stage. At the peak 
of the wave, the relative difference between the Re = 20Af and Re = 100M waveforms is less than 0.6% 
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times when the numerical noise in the waveform is comparable to the true signal amplitude and at very early times 
due to the inability of the codes to accurately model the relatively high-frequency initial data burst. 

Note that the magnitude and phase differences between ip^ extraction and characteristic extraction indicated in 
Figures [T4"! and IT51 are based upon the perturbative extrapolation formula (|5.1[) . It is also common practice to compute 
-04 at large radii and then extrapolate the values to infinity, cf. the waveform comparisons in the report of the Samurai 
project j56|. In carrying out the extrapolation, the waveforms are translated by r* , where r* = r + 2M\og(r/2M — 1) 
is the tortoise coordinate obtained from the areal radius of the extraction sphere r and M is the ADM mass of the 
system 57). Here we use a linear extrapolation based upon waveforms at R = 50M and R — 100M to obtain an 
extrapolation rip^liri) on ,J? + that is accurate to order 0(1/R 2 ). The deviations from the characteristic waveform 
are displayed in Figure[TJJl where we plot Re[5ijji{R — 50M)], Re^ip^R = 100M)] and RelSip^lin)} and in Figure [T71 
where we plot the phase differences 5§(R = 50M), 8$(R = 100M) and 8$(lin). The plots show how the deviations 
decrease with extraction radius. Linear extrapolation considerably reduces the deviation but it is interesting that 
perturbative extrapolation via (I5.1|) . which is based upon the single R — 100 M result, gives the smallest deviation. 



As we discuss next, such time domain comparisons can be of deceptive value for gravitational wave data analysis. 
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FIG. 14: Plots of the time dependence of the real part of the (2,2) spherical harmonic components of 8ip4,{t), as defined in 
(|6.4[) . and the characteristic waveform ^(t). Here 5ip4,(t) measures the difference between the perturbative and characteristic 
values of ^(t). The approximate 1% ratio between the peak amplitudes of 8ip4.(t) and gives an estimate of the difference 
between perturbative and characteristic extraction. 



VII. ADVANCED LIGO ACCURACY STANDARDS 

It has been emphasized [58[ that the direct use of time domain errors would be an abuse of the accuracy standards 
required of model waveforms to be suitable for gravitational wave data analysis. The raw error envelopes SN(t), 5^(t) 
and Sip^t) cannot be used to test whether the accuracy standards are satisfied. Proper accuracy standards must take 
into account the power spectral density of the detector noise S n (f), which is calibrated with respect to the frequency 
domain strain h(f). Consequently, the primary accuracy standards must be formulated in the frequency domain in 
order to take detector sensitivity into account. Fortunately, for the purpose of calibrating waveforms from numerical 
simulations, it has been possible to translate the frequency domain accuracy requirements into requirements on the 
time domain Li error norms which meet all the needed criteria [l], [33] • 

There are two distinct criteria for waveform accuracy. First, if the numerical waveform were not sufficiently accurate 
then an unacceptable fraction of real signals would pass undetected through the corresponding filter. Second, the 
accuracy impacts on whether a detected waveform measures the physical properties of the source, e.g. mass and spin, 
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FIG. 15: Plot of the time dependence of the difference in phase <5$, measured in radians, between the (I = 2, m = 2) components 
of the characteristic waveform ^(t) and the perturbative waveform ripi(t)/2. The phase differences are less than 0.05 radians 
in the interval 40M < t < 250A/ beginning after the initial burst of junk radiation and extending into the late ringdown. The 
phase errors become large at late times when the numerical noise is comparable to the true signal amplitude and at very early 
times due the inability of the codes to accurately model the relatively high-frequency initial data burst. 



to a level commensurate with the accuracy of the observational data. The accuracy standards for model waveforms 
have been formulated to prevent these potential losses in the detection and scientific measurement of gravitational 
waves. 

For a numerical waveform with strain component h(t), the time domain error is measured by 



\\Sh\\ 



(7.1) 



where Sh is the error in the numerical approximation and H-FH 2 = J dt\F(t)\ 2 , i.e. \\F\\ is the L>2 norm, which 
in principle should be integrated over the complete time domain of the model waveform obtained by splicing a 
perturbative chirp waveform to a numerical waveform for the inspiral and merger. 

The error can also be measured in terms of time derivatives of the strain. In our case, the first time derivative 
corresponds to the error in the news 



£i(ReN) 



||&ReiV|| 



S^ImN) 



\\5ImN\\ 



\\ReN\\ ' 1V ' \\ImN\\ 
and the second time derivative corresponds to the Weyl component error 

\SReM „ _ ||<5Im*|| 



£a(ite¥) 



I .Rett I 



(7.2) 



(7.3) 



Here we measure 5N and (5'F according to (|6.2I) and (16.31) and use the 3rd order Richardson extrapolations to compute 
||J?eJV||, ||JmiV||, | jibe's 1 1 and ||/mW||. It is also of interest to measure the "error" 



(7.4) 
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FIG. 16: Plots of the time dependence of the difference Re[5ip4] between the (I = 2, m = 2) components of the characteristic 
waveform ^(t) and the Cauchy ip4 waveforms extracted at R — 50M and R — 100M, and their linear extrapolation to 
R — co (denoted by "lin"). For comparison, we also include the corresponding difference (denoted by "pert") using the 
perturbative extrapolation (|5.1|) . The plots show the expected trend toward smaller errors as R — s- oo. Interestingly, perturbative 
extrapolation, which uses only the R = 100M extraction sphere, gives the smallest deviation. 



corresponding to the difference between perturbative and characteristic extraction, where Sipi is normalized according 

to (ga . 

In [l|, it was shown that sufficient conditions to satisfy data analysis criteria for detection and measurement can 
be formulated in terms of any of the error norms £k = (£o, £i, £2), i-e. in terms of the strain, the news or the Weyl 
component. The accuracy requirement derived in [l| for detection is 



and the requirement for measurement is 



£k < C k \/2e 

1: 



£k < Cfc — . 

p 



Here p is the optimal signal-to-noise ratio of the detector, defined by 

4|M/)| 2 



Sn(f) 



-df; 



(7.5) 



(7.6) 



(7.7) 



Ck are dimensionless factors introduced in [l| to rescale the traditional signal-to-noise ratio p in making the transition 
from frequency domain standards to time domain standards; e ma x determines the fraction of detections lost due 
to template mismatch, cf. Eq. (14) of [37^ ; and rj c < 1 corrects for error introduced in detector calibration. These 
requirements for detection and measurement, for either k = 0, 1, 2, conservatively overstate the basic frequency domain 
requirements by replacing S n (f) by its minimum value in transforming to the time domain. 

The values of Ck for the inspiral and merger of non-spinning, equal-mass black holes have been calculated in [l[ 
for the advanced LIGO noise spectrum. As the total mass of the binary varies from — ¥ 00, Cq varies between 
.65 > Co > 0, C\ varies between .24 < C\ < .8 and C2 varies between < Ci < 1. Thus only the error £1 in the news 
can satisfy the criteria over the entire mass range. The error in the strain £q provides the easiest way to satisfy the 
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FIG. 17: Plots of the time dependence of the phase difference 5$ between the (I = 2, m = 2) components of the characteristic 
waveform ^(t) and the Cauchy tpi waveforms extracted at i? = 50M and i? = 100M, and the corresponding linear extrapolation 
to R = co (denoted by "lin"). For comparison, we also include the corresponding <53> (denoted by "pert") obtained by 
perturbative extrapolation (|5,l|l . The plots show the expected trend toward smaller errors as R — > co. Interestingly, perturbative 
extrapolation, which uses only the R = lOOAf extraction sphere, gives the smallest deviation. 



criteria in the low mass case M « Mq and the error in the Weyl component £2 provides the easiest way to satisfy 
the criteria in the high mass case M >> Mq. 

We first concentrate on the error in the news, for which the accuracy requirement for detection is 



£\ < Ciy/2e max , (7.8) 

and the requirement for measurement is 

£i<Ci— . (7.9) 
P 

Table [V] gives values of several versions of the £ 1 error for the inspiral and merger of non-spinning, equal mass 
black holes described in Sec. |Vl For practical purposes, the error norms were computed over the time period of the 
simulation rather than for a complete model waveform obtained by splicing to a post-Newtonian chirp waveform. 
Assuming that the nonlinear error in the chirp waveform is small compared to the error in the numerical waveform, 
the effect is to overestimate the error norms by underestimating the denominators in (17.11) . (|7.2I) and (|7.3[) . However, 
it has been pointed out in (59l . l60j that splicing to a chirp waveform can produce significant error unless the numerical 
waveform extends to a large number of orbits, which can be computationally prohibitive; otherwise, only for binary 
masses > 100M Q is the splicing error negligible. 

An advantage of the £ 1 error norm based upon the news function is that the denominator in (|T. 1 [) is directly related 
to the radiated energy. As a result, the factor by which £\{N) is overestimated is J 7-1 / 2 , where 

T , = A£(Numerical) 

A£(Chirp) + AS(Numerical) 1 ' ' 



and AE denotes the energy radiated in the indicated time periods. The total energy radiated during the post- 
Newtonian inspiral and merger can be estimated from the difference between the final black hole mass Mh and the 
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mass Mo of the binary for a large initial orbit. The energy Ai^(Numerical) radiated during the numerically modeled 
time period can be obtained from the Bondi mass-loss formula. For the binary inspiral being considered here, the 
final black hole mass is Mh ~ 0.965187; the initial mass of the system at infinite separation (given by the sum of the 
individual black hole masses) is Mq ps 1.01447; and AS (Numerical) « 0.0346. This leads to the fraction of energy 
T s» .702 radiated during the numerical period, or 

T~ 1/2 « 1.19 (7.11) 

for the factor by which the £\ errors in Table [V] are overestimated. We re-emphasize that the values in the Table do 
not include the error introduced by splicing the post-Newtonian and numerical waveforms. 

Besides the values £i(N) of the numerical truncation error in the real and imaginary part of the news function 
extracted at Re = 20M, 50M and 100M, Table IV1 includes the corresponding truncation error fi(iV^) obtained from 
integrating "J via (|3.4p . The Table also includes the modeling errors £i(N^ R ^ 2 o,ioo)) an d ^i(-^a_r(50,ioo)) m the news 
which results from the differences Nr = 2q — -/V#=ioo and N R=50 — -/Vr =10 o obtained from extracting the waveform at 
radii Re = 50M and Re = 20M compared with extraction at Re = 100M. In computing these error norms, we 
integrate over the interval corresponding to t/M > 100 in Fig. [13] to eliminate effects of the initial junk radiation. 



TABLE V: Error norms of the (2, 2) spherical harmonic mode for the Bondi news N, its counterpart (obtained by time 
integral of the Weyl component $) and for the differences Nar comparing extraction at radii Re = 20M and Re = 50M to 
extraction at Re = 100M. 
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Consider first the criterion for detection where we set e max = .005, which for advanced LIGO ensures less than a 
10% signal loss, a target which is often adopted in LIGO searches for compact binaries 37]. For this target, (|7.5[) 
reduces to £\ < O.lCi, or £\ < .024 for the low mass bound C\ ~ .24. This criterion is easily satisfied by all the 
error norms in Table [V] Thus the advanced LIGO detection criterion is satisfied by CCE waveforms obtained from 
either the news or Weyl component throughout the entire binary mass range. In addition, the detection criterion 
is unaffected by modeling errors introduced by choice of extraction radius. Note that the £\{N) and £i(N^) errors 
decrease with larger extraction radius. This is expected since the truncation error introduced by the characteristic 
evolution code depends upon the size of the integration region between the extraction worldtube and . 

The criterion for measurement is more stringent. For a calibration factor given by the expected lower bound 
ijmin = 0.4 and for the lower bound C\ ~ .24 corresponding to the small mass limit, (|7.9j) reduces to 

gi<Ci!j£= 9.6x10-3 
P P 

For the most optimistic advanced LIGO signal-to-noise ratio, which is expected to be p sa 100 for the strongest and 
best tuned events, the requirement for measurement is then £ x < 9.6 x 10~~ 4 . Thus, comparing (|7.12[) to the values 
in Table |V1 the advanced LIGO measurement criterion is satisfied throughout the entire binary mass range by the 
numerical truncation error £i{N) in the CCE waveform obtained from the news function. The £\(N^) error obtained 
from the Weyl component for extraction radii Re > 50M also satisfy this full range of measurement standards. The 
value of £i(Ny) for Re = 20M would satisfy the full range of measurement standards for p < 100 if reduced by the 
factor J 7-1 / 2 given in (|7.11|) . Note also that the truncation error is being conservatively measured by the 0(A 2 ) error 
(|6.3[) , rather than the third order accurate error in the Richardson extrapolated waveform. 

These results can be compared with the measurement criterion for advanced LIGO data analysis reported in the 
Samurai project (5(|, which was also based upon a non-spinning, equal-mass binary black hole inspiral and merger. 
There it was found that the mismatch between perturbative waveforms obtained using various Cauchy codes limited 
the measurement application to signal-to- noise ratios p < 25. This is consistent with our experience, and that reported 
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in that the additional truncation error introduced by applying CCE to a Cauchy simulation of a binary inspiral 
is much smaller than the numerical error resulting from the Cauchy code. 

The values of£i(JV A fl) in Table|V] give an estimate of the modeling error introduced by different choices of extraction 
radius. The error £i(-^Vah(50,ioo))j introduced by extraction at Re = 50M as compared to Re = 100M , only satisfies 
the full range of measurement standards for signal-to-noise ratios p < 21, or p < 25 if (|7.1ip is taken into account. 
This would cover the most likely advanced LIGO events. These modeling errors primarily result from initialization 
effects which we have discussed and which would be less significant in simulations with a higher number of orbits. The 
results suggest that the choice of extraction radius should be balanced between a sufficiently large radius to reduce 
initialization effects and a sufficiently small radius where the Cauchy grid is more highly refined and outer boundary 
effects arc better isolated. For the Cauchy grid setup in the present case, there is a factor of 2 in refinement at 
r = 50M compared to r — 100A/, which for 8th order finite differencing has considerable impact on the error. Future 
experiments with longer runs involving more orbits will supply valuable guidance for optimizing the extraction radius. 

Table PVTl gives some pertinent £2 error norms for the (2,2) spherical harmonic component. Besides the numerical 
truncation error £2^) obtained for characteristic extraction at Re = 20M, Re — 50M and Re = 100M, the Table 
includes the error norm £2(Sip) measuring the difference between perturbative and characteristic extraction, as defined 
in (|7.4p . obtained at Re = 50 Af and Re = 100M. The Table also includes the modeling error f2(^A_R(50,ioo)) resulting 
from the difference ^r = 5o — ^r=ioo obtained using characteristic extraction at radii Re — 50M and Re = 100M, 
as well as the corresponding error norm £2(^4, ah(50, 100)) resulting from the difference obtained using perturbative 
extraction at Re = 50M and Re = 100M, i.e. 

c ( U 1 \ ll^e[(rV> 4 /2)| r= iooM - H> 4 /2)| r = 5 oM]|| 
t2{Reip4,AR{50,ioo)) = ll-Re^ll 

c (T I \ \\ Im l( r ^i/ 2 )\r=W0M ~ M 4 /2)|^50A/]|| ,„ , 

£ 2 (/W> 4 ,A.R(50,100)) = \\Im^\\ ■ ( 7 - 14 ) 

All norms are again computed over the time interval t/M > 100 indicated in Fig. [T3]to reduce effects of initial junk 
radiation. 

Although the £2 norms are not effective for low mass binaries, they give some useful information for comparing 
extraction at various radii and comparing characteristic and perturbative extraction. In the high mass limit for which 
C2 — 1 and with the same lower limits for e m ax and rj c as for the £\ norms, the detection criterion (|7.5|) reduces to 



and the measurement criterion (|7.6p reduces to 



£ 2 <V2e^=0.1 (7.15) 



£2 < V - = ™. (7.16) 
P P 



All the error norms in Table IVII satisfy the detection requirement for this high mass limit. The truncation errors 
£2^) decrease with extraction radius as in the case of the £\{N^,) errors. The values at all three extraction radii 
satisfy the measurement requirement for the most optimistic advanced LIGO signal-to-noise ratio p = 100. These 
results are consistent with the £\{N^) error in Table IVl obtained by integrating 

The norms £2(Sip) measure the difference between characteristic and perturbative extraction. The results in the 
Table show that this difference is fairly independent of whether the waveforms are extracted at Re = 50M or 
Re = lOOAf . In the high mass limit in which (|7.16p is valid, these errors impact the measurement criterion only for 
signal to noise ratios p > 59 but they could be expected to be more significant for low mass binaries. Whether the 
£2(81/}) error can be attributed to characteristic extraction or to perturbative extraction cannot be decided from this 
single test and deserves further investigation. A definitive answer would of course require knowledge of the "exact" 
waveform. 

The modeling error ^(V^Aii^o.ioo))! which results from comparing perturbative extraction at Re = 50M and 
Re = 100M, is considerably larger than the corresponding modeling error £ 2 (^ , Afl(50,ioo)) f° r characteristic extraction. 
This confirms the expectation that perturbative extraction requires a large extraction radius. Both of these modeling 
errors are substantial, which further emphasizes the importance of an optimal choice of extraction radius. 



VIII. CONCLUSION 



We have developed a new characteristic waveform extraction tool. Bugs and inconsistencies in the previous version 
have been eliminated. The extracted waveform from a binary black hole inspiral now shows clean convergence. We have 
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TABLE VI: £2 error norms for the (2, 2) spherical harmonic mode of the CCE waveform ^ obtained using extraction worldtube 
radii Re = 20M, Re = 50 M and Re = 100M and the norms of the difference 5xp between *t and the perturbative ipi waveforms 
extracted at 50M and 100M. We also tabulate the modeling error £2(^^^(50,100)) resulting from the difference in extracting 
\P at 50M and 100M, and the corresponding modeling error £2(V , 4,a«(50,ioo)) f° r extraction via ip4 



Variable 


Re 


Im 


&(*)fl=20 


1.14 x 10" 


3 


1.17 x 10" 


-3 


£ 2 (*)fl=50 


4.04 x 10" 


4 


3.53 x 10" 


-4 


£2(^)11=100 


2.81 x 10" 


4 


2.09 x 10" 


-4 


£2(^)^=50 


5.09 x 10" 


3 


5.08 x 10" 


-3 


£2(5^)fl=100 


6.81 x 10" 


3 


6.32 x 10" 


-3 


£2(*Afl(50,100)) 


1.94 x 10" 


2 


1.91 x 10" 


-2 


£2(^4^^(50,100)) 


3.13 x 10" 


2 


3.14 x 10" 


-2 



demonstrated that this allows the use of Richardson extrapolation to obtain third order accurate waveforms whose 
numerical truncation error satisfies the advanced LIGO standards for detection and measurement. Characteristic 
waveform extraction from a binary black hole inspiral can now be obtained without any recourse to linearization 
and from extraction radii as small as R = 20M. The Cauchy interface has been simplified in terms of a spectral 
decomposition. 

There are still elements where accuracy could be improved. Some of these, such as more accurate start-up algorithms 
for the radial integrations at the extraction worldtube and more accurate asymptotic limits at might be handled 
by small modifications but others, such as extending the overall accuracy to 4th or higher order, would entail a 
more major overhaul of the underlying PITT code. This is perhaps long overdue, but a proper treatment would 
require a better understanding of the underlying mathematical problem. The well-posedness of the gravitational 
worldtube-nullcone initial-boundary value problem upon which the code is based has not yet been established. Only 
recently has well-posedness been demonstrated for the corresponding nonlinear scalar wave problem (6lj . The PITT 
code was developed in the early days of numerical relativity when considerations of well-posedness did not arise in 
the formulation of Cauchy as well as characteristic codes. The development of a stable characteristic code involved 
"educated guesses" . Today, the numerical relativity community is more aware of the benefits that a well-posed problem 
can bring. Most important, a proof of well-posedness of the continuum problem by means of energy estimates can 
be converted to ensure stability of the corresponding finite difference problem by the analogous estimates obtained 
by summation by parts. A new characteristic code based upon this approach would be of great value. Of equal 
value would be the implementation of Cauchy-characteristic matching (CCM), in which the characteristic evolution 
is used to supply outer boundary data for the Cauchy evolution. So far, CCM has only been successfully applied to 
a harmonic Cauchy code in the linearized regime [62| . 

Although there is room for further improvement in the CCE tool presented here, there also is pressing interest 
from several numerical relativity groups to apply the tool to extract waveforms from binary black hole inspirals. The 
emerging importance of this problem to the future of gravitational wave astronomy has created an urgency to make 
characteristic waveform extraction widely available. Simulations of binary black hole inspirals are too computationally 
expensive to be carried out solely for the purpose of wave extraction tests. This would conflict with the demands 
to apply computational resources to results of importance to gravitational wave astronomy and binary black hole 
astrophysics. However, the extra computational expense of adding characteristic extraction to a Cauchy simulation is 
fractionally small. For our tests, where we extracted twice as often as required, the interpolation, decomposition, and 
saving of the metric data used only « 6.9% of the total simulation time. The application of characteristic extraction 
to simulations of astrophysical importance will at the same time provide a practical approach to improving the 
extraction tool by comparing results obtained with different formulations, different numerical techniques and different 
grid specifications. In particular, our test results emphasize the need for a better understanding of the optimal choice 
of extraction radius, which would balance between the discretization error in the Cauchy code, the initialization error, 
the error originating at the outer Cauchy boundary and the relatively small discretization error from CCE. 

We have demonstrated here how the module can be applied to the LazEv code, which is a finite-difference BSSN 
code, to produce calibrated binary black hole waveforms. We welcome applications to codes based upon other 
formulations of the Einstein equations, e.g. the harmonic formulation, and based upon other numerical methods, e.g. 
spectral methods. In particular, characteristic extraction offers a way to unambiguously compare binary black hole 
waveforms obtained from the same initial data using codes based upon different formulations of the Einstein equations, 
different numerical techniques, different evolution gauges and different methods of treating the internal singularities 



29 



(by punctures or by excision). Such comparisons would be of especial importance in the case of precessing binaries 
composed of high spin black holes, where the reliability of perturbative extraction has not been extensively tested. 

We have made the present characteristic waveform extraction tool publicly available as part of the Einstein 
Toolkit [H. 
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Appendix A: Code revision 

Revisions to the worldtube module: 

• The numerical error in the previous version of the worldtube module did not converge properly upon grid 
refinement. We have traced this problem to an inconsistency in the startup algorithm for the integration of 
the characteristic equations away from the extraction worldtube. Data from the Cauchy code had been used 
in an overdetermined manner to supply the integration constants for the characteristic equations. As a result, 
the Cauchy evolution introduced inconsistencies with the characteristic equations which degraded convergence 
of the numerical error. We have revamped this start-up algorithm so that the worldtube module now has clean 
second order accuracy with respect to grid size. 

• We have found and corrected bugs which had been introduced in the implementation of features designed to 
improve code performance. In particular, in parallelizing the code using the Cactus framework [50|, a complex 
spin- weighted term in the evolution module was incorrectly declared to be a real variable. In addition, it had been 
suggested that improved accuracy could be obtained by reducing second derivatives in the angular directions 
to first order form by the use of auxiliary variables [4lj |. In the process of doing so, values of certain variables 
in the subroutines for the data at the extraction worldtube were inadvertently interchanged between the North 
and South stereographic patches. The introduction of these bugs made the resulting code inconsistent with the 
Einstein equations. (From tracing through the code archive, we determined that the bugs were introduced in 
2002 or later so that they do not affect the validity of results prior to 2002. C. Reisswig has informed us that 
he recomputed some of the results in [H, [|| using the corrected code and found good qualitative agreement with 
the original results.) 

• The matching interface has been simplified by introducing a pseudospectral decomposition of the Cauchy metric 
in the neighborhood of the extraction worldtube. This provides more economical storage of the inner boundary 
data for the characteristic code so that the waveform at can be obtained with small computational burden 
compared to the Cauchy evolution. 

• Interpolation error arises because the characteristic grid points do not lie exactly on the extraction worldtube 
determined by the Cauchy coordinates. The interpolation stencils change in a discontinuous way when the grid 
is refined. Consequently, although this interpolation error is second order in grid size, there is a small stochastic 
component relative to the choice of grid. This can obscure the results of convergence tests. We have reduced 
such sources of error so that convergence tests can be used to validate the interface modules. 

• We have streamlined the start-up procedure at the extraction worldtube by initializing the auxiliary variables 
(introduced to remove second angular derivatives) directly in terms of the main variables. 
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• In previous applications of the extraction module, it was expedient to set the characteristic data on the initial 
hypersurface to zero outside some radius. This necessitated a transition region to obtain continuity with the 
initial Cauchy data, which requires non-zero initial characteristic data at the extraction worldtube. Here we 
initialize the data by requiring that the Newman- Penrose component of the Weyl tensor intrinsic to the initial null 
hypersurface vanish, i.e. ipo — 0. For a linear perturbation of the Schwarzschild metric, this condition eliminates 
incoming radiation crossing the initial null hypersurface. Since tpo consists of a second radial derivative of the 
characteristic data, this condition allows both continuity at the extraction worldtube and the desired asymptotic 
falloff of the characteristic data at infinity. 

Modifications of the PITT code: 

• A source of error in characteristic evolution is the intergrid interpolations arising from the stereographic patches 
used to coordinatize the spherical cross-sections of the outgoing null hypersurfaces. The previous version of the 
code used two square stereographic patches centered about the North and South poles, each overlapping the 
equator. This has now been modified by shrinking the overlap region so that each patch has a circular boundary 
located slightly past the equator, as is the practice in the use of stereographic grids in meteorology [64{. This 
eliminates the region near the corners of the square patch where the numerical error was most troublesome. 
Angular numerical dissipation has also been introduced and shown to be effective in controlling the short 
wavelength noise arising from the intergrid interpolations across the stereographic patches. Tests show that the 
resulting waveforms have smooth numerical error as functions on the sphere Q . 



Characteristic codes based upon a six patch covering of the sphere [65, 66] offer the potential for better accuracy 
but they have not yet been developed to handle waveform extraction. See [5j for a comparison of the six patch 
and the stereographic approaches on a test problem. 

• The accuracy of the angular derivatives has been increased to a 4th order finite difference approximation, as 
opposed to the 2nd order accuracy in the original code. The radial derivatives and time integration remain 
second order accurate. 

• Some of the differential equations governing propagation along the characteristics become degenerate at J? + and 
affect the accuracy of asymptotic quantities such as the Bondi news function. The correct asymptotic behavior 
has now been incorporated into the finite difference approximation in order to increase accuracy. In addition, 
the accuracy of certain one-sided finite difference approximations necessary at + has also been improved. 

• In addition, the code has been extended to supply the waveform at in terms of the radiative component of 
the Weyl tensor as well as the Bondi news function. For tests in the linearized regime, extraction via the Weyl 
tensor was found to be slightly more accurate than via the news function when large gauge effects are introduced 
in the characteristic coordinates [5|. On the other hand, the higher derivatives involved in computing the Weyl 
tensor lead to less smoothness in the numerical error. Overall, the two methods are competitive. 
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